function [b,se,score]=PostDouble(y,x,z)

n = numel(y);

% a1 = lasso2step(z,x);
% a2 = lasso2step(z,y);

yt = y - mean(y);
xt = x - mean(x);
zt = z - ones(n,1)*mean(z);

% Takes WAY TOO LONG given that we are looping over each coefficient
% [tmp1,tmp2] = lasso(zt,xt,'CV',5);
% a1 = feasiblePostLasso(xt,zt,'MaxIter',0,'beta0',tmp1(:,tmp2.IndexMinMSE));
% [tmp1,tmp2] = lasso(zt,yt,'CV',5);
% a2 = feasiblePostLasso(yt,zt,'MaxIter',0,'beta0',tmp1(:,tmp2.IndexMinMSE));

% Original
% a1 = feasiblePostLasso(xt,zt,'MaxIter',5);
% a2 = feasiblePostLasso(yt,zt,'MaxIter',5);

% Start from using five most marginally correlated regressors
czx = corr(xt,zt);
czy = corr(yt,zt);
czxs = sort(-abs(czx));
czys = sort(-abs(czy));
suppx0 = find(abs(czx) >= -czxs(5));
suppy0 = find(abs(czy) >= -czys(5));
tmpx0 = zt(:,suppx0)\xt;
tmpy0 = zt(:,suppy0)\yt;
bx0 = zeros(size(zt,2),1);
by0 = zeros(size(zt,2),1);
bx0(suppx0) = tmpx0;
by0(suppy0) = tmpy0;
a1 = feasiblePostLasso(xt,zt,'MaxIter',1,'beta0',bx0);
a2 = feasiblePostLasso(yt,zt,'MaxIter',1,'beta0',by0);


supp1 = find(a1 ~=0);
supp2 = find(a2 ~=0);
supp = union(supp1,supp2);


x1 = [x,z(:,supp),ones(n,1)];
XpXinv = pinv(x1'*x1);
b1 = XpXinv*(x1'*y); 
e1 = y-x1*b1;  
b = b1(1);
se1 = hetero_se(x1,e1,XpXinv);
se = se1(1);

vh = xt - zt(:,supp)*(zt(:,supp)\xt);
uh = yt - zt(:,supp)*(zt(:,supp)\yt);
eh = uh - vh*b;

score = (vh.*eh)/(vh'*vh);

